packages = c('rgdal', 'sf', 'tmap', 'tidyverse', 'sp', 'rgeos','maptools', 'raster', 'spatstat', 'tmaptools', 'spdep', 'OpenStreetMap')
for (p in packages){
if (!require(p, character.only = T)){
install.packages(p)
}
library(p, character.only = T)
}
## Loading required package: rgdal
## Loading required package: sp
## rgdal: version: 1.5-16, (SVN revision 1050)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 3.0.4, released 2020/01/28
## Path to GDAL shared files: C:/Users/keith/Documents/R/win-library/4.0/rgdal/gdal
## GDAL binary built with GEOS: TRUE
## Loaded PROJ runtime: Rel. 6.3.1, February 10th, 2020, [PJ_VERSION: 631]
## Path to PROJ shared files: C:/Users/keith/Documents/R/win-library/4.0/rgdal/proj
## Linking to sp version:1.4-2
## To mute warnings of possible GDAL/OSR exportToProj4() degradation,
## use options("rgdal_show_exportToProj4_warnings"="none") before loading rgdal.
## Loading required package: sf
## Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 6.3.1
## Loading required package: tmap
## Loading required package: tidyverse
## -- Attaching packages -------------------------------------- tidyverse 1.3.0 --
## v ggplot2 3.3.2 v purrr 0.3.4
## v tibble 3.0.3 v dplyr 1.0.1
## v tidyr 1.1.1 v stringr 1.4.0
## v readr 1.3.1 v forcats 0.5.0
## -- Conflicts ----------------------------------------- tidyverse_conflicts() --
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()
## Loading required package: rgeos
## rgeos version: 0.5-3, (SVN revision 634)
## GEOS runtime version: 3.8.0-CAPI-1.13.1
## Linking to sp version: 1.4-2
## Polygon checking: TRUE
## Loading required package: maptools
## Checking rgeos availability: TRUE
## Loading required package: raster
##
## Attaching package: 'raster'
## The following object is masked from 'package:dplyr':
##
## select
## The following object is masked from 'package:tidyr':
##
## extract
## Loading required package: spatstat
## Loading required package: spatstat.data
## Loading required package: nlme
##
## Attaching package: 'nlme'
## The following object is masked from 'package:raster':
##
## getData
## The following object is masked from 'package:dplyr':
##
## collapse
## Loading required package: rpart
## Registered S3 method overwritten by 'spatstat':
## method from
## print.boxx cli
##
## spatstat 1.64-1 (nickname: 'Help you I can, yes!')
## For an introduction to spatstat, type 'beginner'
##
## Note: spatstat version 1.64-1 is out of date by more than 5 months; we recommend upgrading to the latest version.
##
## Attaching package: 'spatstat'
## The following objects are masked from 'package:raster':
##
## area, rotate, shift
## Loading required package: tmaptools
## Loading required package: spdep
## Loading required package: spData
## To access larger datasets in this package, install the spDataLarge
## package with: `install.packages('spDataLarge',
## repos='https://nowosad.github.io/drat/', type='source')`
## Loading required package: OpenStreetMap
popdata <- read_csv('data/planning-area-subzone-age-group-sex-and-type-of-dwelling-june-2011-2019.csv')
## Parsed with column specification:
## cols(
## planning_area = col_character(),
## subzone = col_character(),
## age_group = col_character(),
## sex = col_character(),
## type_of_dwelling = col_character(),
## resident_count = col_double(),
## year = col_double()
## )
popdata2019 <- popdata %>%
filter(year == 2019) %>%
filter(age_group == "65_to_69" | age_group == "70_to_74" | age_group == '75_to_79' | age_group == '80_to_84' | age_group == '85_to_89' | age_group == '90_and_over') %>%
group_by(planning_area, subzone, age_group) %>%
summarise(elderly_count = sum(resident_count)) %>%
ungroup() %>%
spread(age_group, elderly_count) %>%
mutate(`elderly_count` = `65_to_69` + `70_to_74` + `75_to_79` + `80_to_84` + `85_to_89` + `90_and_over`)
## `summarise()` regrouping output by 'planning_area', 'subzone' (override with `.groups` argument)
popdata2019 <- mutate_at(popdata2019, .vars = c("subzone", "planning_area"), .funs=toupper)
eldercare_sf <- st_read(dsn='data', layer='ELDERCARE')
## Reading layer `ELDERCARE' from data source `C:\Users\keith\Desktop\IS415 Geospatial Analytics\SPP\SPP\data' using driver `ESRI Shapefile'
## Simple feature collection with 133 features and 18 fields
## geometry type: POINT
## dimension: XY
## bbox: xmin: 14481.92 ymin: 28218.43 xmax: 41665.14 ymax: 46804.9
## projected CRS: SVY21
eldercare_sf <- eldercare_sf %>%
mutate(label = "Eldercare centres")
Assign EPSG
eldercare_sf <- st_transform(eldercare_sf, 3414)
st_crs(eldercare_sf)
## Coordinate Reference System:
## User input: EPSG:3414
## wkt:
## PROJCRS["SVY21 / Singapore TM",
## BASEGEOGCRS["SVY21",
## DATUM["SVY21",
## ELLIPSOID["WGS 84",6378137,298.257223563,
## LENGTHUNIT["metre",1]]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433]],
## ID["EPSG",4757]],
## CONVERSION["Singapore Transverse Mercator",
## METHOD["Transverse Mercator",
## ID["EPSG",9807]],
## PARAMETER["Latitude of natural origin",1.36666666666667,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8801]],
## PARAMETER["Longitude of natural origin",103.833333333333,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8802]],
## PARAMETER["Scale factor at natural origin",1,
## SCALEUNIT["unity",1],
## ID["EPSG",8805]],
## PARAMETER["False easting",28001.642,
## LENGTHUNIT["metre",1],
## ID["EPSG",8806]],
## PARAMETER["False northing",38744.572,
## LENGTHUNIT["metre",1],
## ID["EPSG",8807]]],
## CS[Cartesian,2],
## AXIS["northing (N)",north,
## ORDER[1],
## LENGTHUNIT["metre",1]],
## AXIS["easting (E)",east,
## ORDER[2],
## LENGTHUNIT["metre",1]],
## USAGE[
## SCOPE["unknown"],
## AREA["Singapore"],
## BBOX[1.13,103.59,1.47,104.07]],
## ID["EPSG",3414]]
Convert sf to sp
eldercare_sp <- as_Spatial(eldercare_sf)
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO"): Discarded datum Unknown based on WGS84 ellipsoid in CRS definition,
## but +towgs84= values preserved
eldercare_spatialpoint <- as(eldercare_sp, "SpatialPoints")
infocomm_sf <- st_read(dsn='data', layer='SILVERINFOCOMM')
## Reading layer `SILVERINFOCOMM' from data source `C:\Users\keith\Desktop\IS415 Geospatial Analytics\SPP\SPP\data' using driver `ESRI Shapefile'
## Simple feature collection with 101 features and 18 fields
## geometry type: POINT
## dimension: XY
## bbox: xmin: 12024.87 ymin: 28385.78 xmax: 41300.53 ymax: 47763.05
## projected CRS: SVY21
infocomm_sf <- infocomm_sf %>%
mutate(label = "Silver Infocomm Junc")
Assign EPSG
infocomm_sf <- st_transform(infocomm_sf, 3414)
st_crs(infocomm_sf)
## Coordinate Reference System:
## User input: EPSG:3414
## wkt:
## PROJCRS["SVY21 / Singapore TM",
## BASEGEOGCRS["SVY21",
## DATUM["SVY21",
## ELLIPSOID["WGS 84",6378137,298.257223563,
## LENGTHUNIT["metre",1]]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433]],
## ID["EPSG",4757]],
## CONVERSION["Singapore Transverse Mercator",
## METHOD["Transverse Mercator",
## ID["EPSG",9807]],
## PARAMETER["Latitude of natural origin",1.36666666666667,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8801]],
## PARAMETER["Longitude of natural origin",103.833333333333,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8802]],
## PARAMETER["Scale factor at natural origin",1,
## SCALEUNIT["unity",1],
## ID["EPSG",8805]],
## PARAMETER["False easting",28001.642,
## LENGTHUNIT["metre",1],
## ID["EPSG",8806]],
## PARAMETER["False northing",38744.572,
## LENGTHUNIT["metre",1],
## ID["EPSG",8807]]],
## CS[Cartesian,2],
## AXIS["northing (N)",north,
## ORDER[1],
## LENGTHUNIT["metre",1]],
## AXIS["easting (E)",east,
## ORDER[2],
## LENGTHUNIT["metre",1]],
## USAGE[
## SCOPE["unknown"],
## AREA["Singapore"],
## BBOX[1.13,103.59,1.47,104.07]],
## ID["EPSG",3414]]
Convert sf to sp
infocomm_sp <- as_Spatial(infocomm_sf)
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO"): Discarded datum Unknown based on WGS84 ellipsoid in CRS definition,
## but +towgs84= values preserved
infocomm_spatialpoint <- as(infocomm_sp, "SpatialPoints")
chas_sf <- st_read(dsn='data/chas-clinics-kml.kml')
## Reading layer `MOH_CHAS_CLINICS' from data source `C:\Users\keith\Desktop\IS415 Geospatial Analytics\SPP\SPP\data\chas-clinics-kml.kml' using driver `KML'
## Simple feature collection with 1167 features and 2 fields
## geometry type: POINT
## dimension: XYZ
## bbox: xmin: 103.5818 ymin: 1.016264 xmax: 103.9903 ymax: 1.456037
## z_range: zmin: 0 zmax: 0
## geographic CRS: WGS 84
chas_sf <- chas_sf %>%
mutate(label = "Chas Clinics") %>%
mutate(capacity = 1)
Assign EPSG
chas_sf <- st_transform(chas_sf, 3414)
st_crs(chas_sf)
## Coordinate Reference System:
## User input: EPSG:3414
## wkt:
## PROJCRS["SVY21 / Singapore TM",
## BASEGEOGCRS["SVY21",
## DATUM["SVY21",
## ELLIPSOID["WGS 84",6378137,298.257223563,
## LENGTHUNIT["metre",1]]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433]],
## ID["EPSG",4757]],
## CONVERSION["Singapore Transverse Mercator",
## METHOD["Transverse Mercator",
## ID["EPSG",9807]],
## PARAMETER["Latitude of natural origin",1.36666666666667,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8801]],
## PARAMETER["Longitude of natural origin",103.833333333333,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8802]],
## PARAMETER["Scale factor at natural origin",1,
## SCALEUNIT["unity",1],
## ID["EPSG",8805]],
## PARAMETER["False easting",28001.642,
## LENGTHUNIT["metre",1],
## ID["EPSG",8806]],
## PARAMETER["False northing",38744.572,
## LENGTHUNIT["metre",1],
## ID["EPSG",8807]]],
## CS[Cartesian,2],
## AXIS["northing (N)",north,
## ORDER[1],
## LENGTHUNIT["metre",1]],
## AXIS["easting (E)",east,
## ORDER[2],
## LENGTHUNIT["metre",1]],
## USAGE[
## SCOPE["unknown"],
## AREA["Singapore"],
## BBOX[1.13,103.59,1.47,104.07]],
## ID["EPSG",3414]]
Convert sf to sp
chas_sp <- as_Spatial(chas_sf)
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO"): Discarded datum Unknown based on WGS84 ellipsoid in CRS definition,
## but +towgs84= values preserved
chas_spatialpoint <- as(chas_sp, "SpatialPoints")
mpsz_sf <- st_read(dsn='data', layer='MP14_SUBZONE_WEB_PL')
## Reading layer `MP14_SUBZONE_WEB_PL' from data source `C:\Users\keith\Desktop\IS415 Geospatial Analytics\SPP\SPP\data' using driver `ESRI Shapefile'
## Simple feature collection with 323 features and 15 fields
## geometry type: MULTIPOLYGON
## dimension: XY
## bbox: xmin: 2667.538 ymin: 15748.72 xmax: 56396.44 ymax: 50256.33
## projected CRS: SVY21
mpsz_sf <- mpsz_sf %>%
dplyr::select(SUBZONE_N, PLN_AREA_N, REGION_N, X_ADDR, Y_ADDR, SHAPE_Leng, SHAPE_Area, geometry)
Assign EPSG
mpsz_sf <- st_transform(mpsz_sf, 3414)
st_crs(mpsz_sf)
## Coordinate Reference System:
## User input: EPSG:3414
## wkt:
## PROJCRS["SVY21 / Singapore TM",
## BASEGEOGCRS["SVY21",
## DATUM["SVY21",
## ELLIPSOID["WGS 84",6378137,298.257223563,
## LENGTHUNIT["metre",1]]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433]],
## ID["EPSG",4757]],
## CONVERSION["Singapore Transverse Mercator",
## METHOD["Transverse Mercator",
## ID["EPSG",9807]],
## PARAMETER["Latitude of natural origin",1.36666666666667,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8801]],
## PARAMETER["Longitude of natural origin",103.833333333333,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8802]],
## PARAMETER["Scale factor at natural origin",1,
## SCALEUNIT["unity",1],
## ID["EPSG",8805]],
## PARAMETER["False easting",28001.642,
## LENGTHUNIT["metre",1],
## ID["EPSG",8806]],
## PARAMETER["False northing",38744.572,
## LENGTHUNIT["metre",1],
## ID["EPSG",8807]]],
## CS[Cartesian,2],
## AXIS["northing (N)",north,
## ORDER[1],
## LENGTHUNIT["metre",1]],
## AXIS["easting (E)",east,
## ORDER[2],
## LENGTHUNIT["metre",1]],
## USAGE[
## SCOPE["unknown"],
## AREA["Singapore"],
## BBOX[1.13,103.59,1.47,104.07]],
## ID["EPSG",3414]]
centroids <- st_centroid(mpsz_sf)
## Warning in st_centroid.sf(mpsz_sf): st_centroid assumes attributes are constant
## over geometries of x
sg <- readOGR(dsn = "data", layer="CostalOutline")
## OGR data source with driver: ESRI Shapefile
## Source: "C:\Users\keith\Desktop\IS415 Geospatial Analytics\SPP\SPP\data", layer: "CostalOutline"
## with 60 features
## It has 4 fields
sg_spatialpoint <- as(sg, "SpatialPolygons")
Joining popdata2019 and mpsz_sf
mpsz_demand <- left_join(mpsz_sf, popdata2019, by=c('PLN_AREA_N' = 'planning_area', 'SUBZONE_N' = 'subzone'))
summary(mpsz_demand)
## SUBZONE_N PLN_AREA_N REGION_N X_ADDR
## Length:323 Length:323 Length:323 Min. : 5093
## Class :character Class :character Class :character 1st Qu.:21864
## Mode :character Mode :character Mode :character Median :28465
## Mean :27257
## 3rd Qu.:31674
## Max. :50425
## Y_ADDR SHAPE_Leng SHAPE_Area 65_to_69
## Min. :19579 Min. : 871.5 Min. : 39438 Min. : 0.0
## 1st Qu.:31776 1st Qu.: 3709.6 1st Qu.: 628261 1st Qu.: 0.0
## Median :35113 Median : 5211.9 Median : 1229894 Median : 270.0
## Mean :36106 Mean : 6524.4 Mean : 2420882 Mean : 686.7
## 3rd Qu.:39869 3rd Qu.: 6942.6 3rd Qu.: 2106483 3rd Qu.:1060.0
## Max. :49553 Max. :68083.9 Max. :69748299 Max. :8120.0
## 70_to_74 75_to_79 80_to_84 85_to_89
## Min. : 0.0 Min. : 0.0 Min. : 0.0 Min. : 0
## 1st Qu.: 0.0 1st Qu.: 0.0 1st Qu.: 0.0 1st Qu.: 0
## Median : 180.0 Median : 110.0 Median : 70.0 Median : 40
## Mean : 466.4 Mean : 294.6 Mean : 193.4 Mean :105
## 3rd Qu.: 760.0 3rd Qu.: 485.0 3rd Qu.: 310.0 3rd Qu.:170
## Max. :4800.0 Max. :2680.0 Max. :1760.0 Max. :970
## 90_and_over elderly_count geometry
## Min. : 0.00 Min. : 0 MULTIPOLYGON :323
## 1st Qu.: 0.00 1st Qu.: 0 epsg:3414 : 0
## Median : 30.00 Median : 730 +proj=tmer...: 0
## Mean : 59.81 Mean : 1806
## 3rd Qu.: 90.00 3rd Qu.: 3000
## Max. :530.00 Max. :18860
mpsz_demand_sp <- as_Spatial(mpsz_demand)
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO"): Discarded datum Unknown based on WGS84 ellipsoid in CRS definition,
## but +towgs84= values preserved
mpsz_demand_spatialpoint <- as(mpsz_demand_sp, "SpatialPolygons")
Calculating Chas Clinics count
mpsz_demand$`chas_count` <- lengths(st_intersects(mpsz_sf,chas_sf))
Calculating eldercare count
mpsz_demand$`eldercare_count` <- lengths(st_intersects(mpsz_sf,eldercare_sf))
Calculating infocomm count
mpsz_demand$`infocomm_count` <- lengths(st_intersects(mpsz_sf,infocomm_sf))
Convert sf to sp
mpsz_sp <- as_Spatial(mpsz_sf)
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO"): Discarded datum Unknown based on WGS84 ellipsoid in CRS definition,
## but +towgs84= values preserved
mpsz_spatialpoint <- as(mpsz_sp, "SpatialPolygons")
sg_osm <- read_osm(mpsz_spatialpoint, ext=1.3)
Class interval for boxmap
boxbreaks <- function(v,mult=1.5) {
qv <- unname(quantile(v))
iqr <- qv[4] - qv[2]
upfence <- qv[4] + mult * iqr
lofence <- qv[2] - mult * iqr
# initialize break points vector. might need more length for bigger data
bb <- vector(mode="numeric",length=7)
# logic for lower and upper fences
if (lofence < qv[1]) { # no lower outliers
bb[1] <- lofence
bb[2] <- floor(qv[1])
} else {
bb[2] <- lofence
bb[1] <- qv[1]
}
if (upfence > qv[5]) { # no upper outliers
bb[7] <- upfence
bb[6] <- ceiling(qv[5])
} else {
bb[6] <- upfence
bb[7] <- qv[5]
}
bb[3:5] <- qv[2:4]
return(bb)
}
get.var <- function(vname,df) {
v <- df[vname] %>% st_set_geometry(NULL)
v <- unname(v[,1])
return(v)
}
boxmap <- function(vnam,df,legtitle=NA,mtitle="Box Map",mult=1.5){
var <- get.var(vnam,df)
bb <- boxbreaks(var)
tm_shape(df) +
tm_fill(vnam,title=legtitle,breaks=bb,palette="-RdBu",
labels = c("lower outlier", "< 25%", "25% - 50%", "50% - 75%","> 75%", "upper outlier")) +
tm_borders() +
tm_layout(title = mtitle, title.position = c("right","bottom"))
}
boxmap("elderly_count", mpsz_demand, mtitle="Elderly Distribution")
Demand and supply of eldercare facilities
tmap_mode("view")
## tmap mode set to interactive viewing
tm_shape(mpsz_demand) +
tm_fill(c("elderly_count"),
title = "Elderly Population",
breaks = c(0, 2500 ,5000, 7500, 10000, 12500, 15000, 17500, 20000),
palette = "Blues",
popup.vars=c("Subzone Name"="SUBZONE_N",
"Planning Area Name"="PLN_AREA_N",
"Elderly Population"="elderly_count"),
showNA = FALSE) +
tm_borders(lwd = 0.1, alpha = 1) +
tm_shape(chas_sf) + tm_dots("red") +
tm_shape(eldercare_sf) + tm_dots("blue")+
tm_shape(infocomm_sf) + tm_dots("green")
Extracting study areas(Tampines, Ang Mo Kio, Serangoon, Jurong West) in spatialpoint dataframe
tp <- mpsz_sp[mpsz_sp@data$PLN_AREA_N == "TAMPINES",]
amk <- mpsz_sp[mpsz_sp@data$PLN_AREA_N == "ANG MO KIO",]
bs <- mpsz_sp[mpsz_sp@data$PLN_AREA_N == "BISHAN",]
jw <- mpsz_sp[mpsz_sp@data$PLN_AREA_N == "JURONG WEST",]
Convert them into generic spatialpolygons objects
tp_spatialpoint <- as(tp, "SpatialPolygons")
amk_spatialpoint <- as(amk, "SpatialPolygons")
bs_spatialpoint <- as(bs, "SpatialPolygons")
jw_spatialpoint <- as(jw, "SpatialPolygons")
creating own objects based on subzones
tp_owin <- as(tp, "owin")
amk_owin <- as(amk, "owin")
bs_owin <- as(bs, "owin")
jw_owin <- as(jw, "owin")
Converting spatialpoints to spatstat ppp
chas_ppp <- as(chas_spatialpoint, "ppp")
infocomm_ppp <- as(infocomm_spatialpoint, "ppp")
eldercare_ppp <- as(eldercare_spatialpoint, "ppp")
using jittering to better visualize overlapped data
chas_ppp_jit <- rjitter(chas_ppp, retry=TRUE, nsim=1, drop=TRUE)
infocomm_ppp_jit <- rjitter(infocomm_ppp, retry=TRUE, nsim=1, drop=TRUE)
eldercare_ppp_jit <- rjitter(eldercare_ppp, retry=TRUE, nsim=1, drop=TRUE)
Combining Chas, Infocomm and Eldercare
ppp_chas_tp <- chas_ppp_jit[tp_owin]
ppp_chas_amk <- chas_ppp_jit[amk_owin]
ppp_chas_bs <- chas_ppp_jit[bs_owin]
ppp_chas_jw <- chas_ppp_jit[jw_owin]
ppp_infocomm_tp <- infocomm_ppp_jit[tp_owin]
ppp_infocomm_amk <- infocomm_ppp_jit[amk_owin]
ppp_infocomm_bs <- infocomm_ppp_jit[bs_owin]
ppp_infocomm_jw <- infocomm_ppp_jit[jw_owin]
ppp_eldercare_tp <- eldercare_ppp_jit[tp_owin]
ppp_eldercare_amk <- eldercare_ppp_jit[amk_owin]
ppp_eldercare_bs <- eldercare_ppp_jit[bs_owin]
ppp_eldercare_jw <- eldercare_ppp_jit[jw_owin]
Convert data to KM by rescaling
ppp_chas_tp.km = rescale(ppp_chas_tp, 1000, "km")
ppp_chas_amk.km = rescale(ppp_chas_amk, 1000, "km")
ppp_chas_bs.km = rescale(ppp_chas_bs, 1000, "km")
ppp_chas_jw.km = rescale(ppp_chas_jw, 1000, "km")
ppp_infocomm_tp.km = rescale(ppp_infocomm_tp, 1000, "km")
ppp_infocomm_amk.km = rescale(ppp_infocomm_amk, 1000, "km")
ppp_infocomm_bs.km = rescale(ppp_infocomm_bs, 1000, "km")
ppp_infocomm_jw.km = rescale(ppp_infocomm_jw, 1000, "km")
ppp_eldercare_tp.km = rescale(ppp_eldercare_tp, 1000, "km")
ppp_eldercare_amk.km = rescale(ppp_eldercare_amk, 1000, "km")
ppp_eldercare_bs.km = rescale(ppp_eldercare_bs, 1000, "km")
ppp_eldercare_jw.km = rescale(ppp_eldercare_jw, 1000, "km")
create owin objects that is required by spatstat.
sg_owin <- as(mpsz_sp, "owin")
sg_owin
## window: polygonal boundary
## enclosing rectangle: [2667.54, 56396.44] x [15748.72, 50256.33] units
Extracting all facilities within the boundaries of the study area
ppp_chas_mpsz = chas_ppp_jit[sg_owin]
ppp_infocomm_mpsz = infocomm_ppp_jit[sg_owin]
ppp_eldercare_mpsz = eldercare_ppp_jit[sg_owin]
Rescaling in KM
ppp_chas_mpsz.km = rescale(ppp_chas_mpsz, 1000, "km")
ppp_infocomm_mpsz.km = rescale(ppp_infocomm_mpsz, 1000, "km")
ppp_eldercare_mpsz.km = rescale(ppp_eldercare_mpsz, 1000, "km")
In order to confirm the spatial patterns for chas clinics, we will have to conduct a hypothesis test for all subzones and chas clinic locations.
H0 = The distribution of chas clinics are randomly distributed. H1= The distribution of chas clinics are not randomly distributed.
The null hypothesis will be rejected if p-value is smaller than alpha < 0.001i
converting spatialdata into ppp object format
performing CSR test
K_chas_tp = Lest(ppp_chas_tp, correction = "Ripley")
K_chas_amk = Lest(ppp_chas_amk, correction = "Ripley")
K_chas_bs = Lest(ppp_chas_bs, correction = "Ripley")
K_chas_jw = Lest(ppp_chas_jw, correction = "Ripley")
K_infocomm_tp = Lest(ppp_infocomm_tp, correction = "Ripley")
K_infocomm_amk = Lest(ppp_infocomm_amk, correction = "Ripley")
K_infocomm_bs = Lest(ppp_infocomm_bs, correction = "Ripley")
K_infocomm_jw = Lest(ppp_infocomm_jw, correction = "Ripley")
K_eldercare_tp = Lest(ppp_eldercare_tp, correction = "Ripley")
K_eldercare_amk = Lest(ppp_eldercare_tp, correction = "Ripley")
K_eldercare_bs = Lest(ppp_eldercare_tp, correction = "Ripley")
K_eldercare_jw = Lest(ppp_eldercare_tp, correction = "Ripley")
K_chas_tp.csr <- envelope(ppp_chas_tp, Kest, nsim = 99, rank = 1, glocal=TRUE)
## Generating 99 simulations of CSR ...
## 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40,
## 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80,
## 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 99.
##
## Done.
K_chas_amk.csr <- envelope(ppp_chas_amk, Kest, nsim = 99, rank = 1, glocal=TRUE)
## Generating 99 simulations of CSR ...
## 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40,
## 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80,
## 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 99.
##
## Done.
K_chas_bs.csr <- envelope(ppp_chas_bs, Kest, nsim = 99, rank = 1, glocal=TRUE)
## Generating 99 simulations of CSR ...
## 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40,
## 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80,
## 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 99.
##
## Done.
K_chas_jw.csr <- envelope(ppp_chas_jw, Kest, nsim = 99, rank = 1, glocal=TRUE)
## Generating 99 simulations of CSR ...
## 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40,
## 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80,
## 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 99.
##
## Done.
K_infocomm_tp.csr <- envelope(ppp_infocomm_tp, Kest, nsim = 99, rank = 1, glocal=TRUE)
## Generating 99 simulations of CSR ...
## 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40,
## 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80,
## 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 99.
##
## Done.
K_infocomm_amk.csr <- envelope(ppp_infocomm_amk, Kest, nsim = 99, rank = 1, glocal=TRUE)
## Generating 99 simulations of CSR ...
## 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40,
## 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80,
## 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 99.
##
## Done.
K_infocomm_bs.csr <- envelope(ppp_infocomm_bs, Kest, nsim = 99, rank = 1, glocal=TRUE)
## Generating 99 simulations of CSR ...
## 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40,
## 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80,
## 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 99.
##
## Done.
K_infocomm_jw.csr <- envelope(ppp_infocomm_jw, Kest, nsim = 99, rank = 1, glocal=TRUE)
## Generating 99 simulations of CSR ...
## 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40,
## 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80,
## 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 99.
##
## Done.
K_eldercare_tp.csr <- envelope(ppp_eldercare_tp, Kest, nsim = 99, rank = 1, glocal=TRUE)
## Generating 99 simulations of CSR ...
## 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40,
## 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80,
## 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 99.
##
## Done.
K_eldercare_amk.csr <- envelope(ppp_eldercare_tp, Kest, nsim = 99, rank = 1, glocal=TRUE)
## Generating 99 simulations of CSR ...
## 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40,
## 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80,
## 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 99.
##
## Done.
K_eldercare_bs.csr <- envelope(ppp_eldercare_tp, Kest, nsim = 99, rank = 1, glocal=TRUE)
## Generating 99 simulations of CSR ...
## 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40,
## 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80,
## 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 99.
##
## Done.
K_eldercare_jw.csr <- envelope(ppp_eldercare_tp, Kest, nsim = 99, rank = 1, glocal=TRUE)
## Generating 99 simulations of CSR ...
## 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40,
## 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80,
## 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 99.
##
## Done.
par(mfrow=c(4,3))
par(mar = rep(3, 4))
plot(K_chas_tp, . - r ~ r,
xlab="d", ylab="L(d)-r",
main="Chas Clinics in Tampines")
plot(K_infocomm_tp.csr, . - r ~ r,
xlab="d", ylab="L(d)-r",
main="Silver infocomm centres in Tampines")
plot(K_eldercare_tp.csr, . - r ~ r,
xlab="d", ylab="L(d)-r",
main="Eldercare services in Tampines")
plot(K_chas_amk.csr, . - r ~ r,
xlab="d", ylab="L(d)-r",
main="Chas Clinics in Ang Mo Kio")
plot(K_infocomm_amk.csr, . - r ~ r,
xlab="d", ylab="L(d)-r",
main="Silver infocomm centres in Ang Mo Kio")
plot(K_eldercare_amk.csr, . - r ~ r,
xlab="d", ylab="L(d)-r",
main="Eldercare services in Ang Mo Kio")
plot(K_chas_bs.csr, . - r ~ r,
xlab="d", ylab="L(d)-r",
main="Chas Clinics in Bishan")
plot(K_infocomm_bs.csr, . - r ~ r,
xlab="d", ylab="L(d)-r",
main="Silver infocomm centres in Bishan")
plot(K_eldercare_bs.csr, . - r ~ r,
xlab="d", ylab="L(d)-r",
main="Eldercare services in Bishan")
plot(K_chas_jw.csr, . - r ~ r,
xlab="d", ylab="L(d)-r",
main="Chas Clinics in Jurong West")
plot(K_infocomm_jw.csr, . - r ~ r,
xlab="d", ylab="L(d)-r",
main="Silver infocomm centres in Jurong West")
plot(K_eldercare_jw.csr, . - r ~ r,
xlab="d", ylab="L(d)-r",
main="Eldercare services in Jurong West")
##Kernal density Estimation with OSM
Create a binding box for selected subzones
tp_bb <- st_bbox(mpsz_sf %>% filter(PLN_AREA_N == "TAMPINES"))
amk_bb <- st_bbox(mpsz_sf %>% filter(PLN_AREA_N == "ANG MO KIO"))
bs_bb <- st_bbox(mpsz_sf %>% filter(PLN_AREA_N == "BISHAN"))
jw_bb <- st_bbox(mpsz_sf %>% filter(PLN_AREA_N == "JURONG WEST"))
Getting osm for the indivisual planning areas
tp_osm <- read_osm(tp_bb, ext=1.1)
amk_osm <- read_osm(amk_bb, ext=1.1)
bs_osm <- read_osm(bs_bb, ext=1.1)
jw_osm <- read_osm(jw_bb, ext=1.1)
getPlnAreaKernelDensityMap <- function(osm, pln, ppp, ppp_str) {
ppp.km <- rescale(ppp, 1000, "km")
kde <- density(ppp.km, sigma=0.20, edge=TRUE, kernel="gaussian")
gridded_kde_bw <- as.SpatialGridDataFrame.im(kde)
kde_bw_raster <- raster(gridded_kde_bw)
projection(kde_bw_raster) <- crs("+init=EPSG:3414 +datum=WGS84 +units=km")
# Plot kernel density map on openstreetmap
tm_shape(osm)+
tm_layout(legend.outside = TRUE, title=ppp_str)+
tm_rgb()+
tm_shape(pln)+
tm_borders(col = "darkblue", lwd = 2, lty="longdash")+
tm_shape(kde_bw_raster) +
tm_raster("v", alpha=0.5,
palette = "YlOrRd")
}
kde_chas_tp <- getPlnAreaKernelDensityMap(tp_osm, tp, ppp_chas_tp, "Chas Clinics in Tampines")
kde_infocomm_tp <- getPlnAreaKernelDensityMap(tp_osm, tp, ppp_infocomm_tp, "Infocomm Centres in Tampines")
kde_eldercare_tp <- getPlnAreaKernelDensityMap(tp_osm, tp, ppp_eldercare_tp, "Eldercare Services in Tampines")
kde_chas_amk <- getPlnAreaKernelDensityMap(amk_osm, amk, ppp_chas_amk, "Chas Clinics in Ang Mo Kio")
kde_infocomm_amk <- getPlnAreaKernelDensityMap(amk_osm, amk, ppp_infocomm_amk, "Infocomm Centres in Ang Mo Kio")
kde_eldercare_amk <- getPlnAreaKernelDensityMap(amk_osm, amk, ppp_eldercare_amk, "Eldercare Services in Ang Mo Kio")
kde_chas_bs <- getPlnAreaKernelDensityMap(bs_osm, bs, ppp_chas_bs, "Chas Clinics in Bishan")
kde_infocomm_bs <- getPlnAreaKernelDensityMap(bs_osm, bs, ppp_infocomm_bs, "Infocomm Centres in Bishan")
kde_eldercare_bs <- getPlnAreaKernelDensityMap(bs_osm, bs, ppp_eldercare_bs, "Eldercare Services in Bishan")
kde_chas_jw <- getPlnAreaKernelDensityMap(jw_osm, jw, ppp_chas_jw, "Chas Clinics in Jurong West")
kde_infocomm_jw <- getPlnAreaKernelDensityMap(jw_osm, jw, ppp_infocomm_jw, "Infocomm Centres in Jurong West")
kde_eldercare_jw <- getPlnAreaKernelDensityMap(jw_osm, jw, ppp_eldercare_jw, "Eldercare Services in Jurong West")
tmap_arrange(kde_chas_tp, kde_infocomm_tp, kde_eldercare_tp, asp=1, ncol=3)
## OpenStreetMapData read with read_osm is static, so not usable in view mode. Please use tm_basemap or tm_tiles, with the provider name set to "OpenStreetMap.Mapnik"
## OpenStreetMapData read with read_osm is static, so not usable in view mode. Please use tm_basemap or tm_tiles, with the provider name set to "OpenStreetMap.Mapnik"
## OpenStreetMapData read with read_osm is static, so not usable in view mode. Please use tm_basemap or tm_tiles, with the provider name set to "OpenStreetMap.Mapnik"
tmap_arrange(kde_chas_amk, kde_infocomm_amk, kde_eldercare_amk, asp=1, ncol=3)
## OpenStreetMapData read with read_osm is static, so not usable in view mode. Please use tm_basemap or tm_tiles, with the provider name set to "OpenStreetMap.Mapnik"
## OpenStreetMapData read with read_osm is static, so not usable in view mode. Please use tm_basemap or tm_tiles, with the provider name set to "OpenStreetMap.Mapnik"
## OpenStreetMapData read with read_osm is static, so not usable in view mode. Please use tm_basemap or tm_tiles, with the provider name set to "OpenStreetMap.Mapnik"
tmap_arrange(kde_chas_bs, kde_infocomm_bs, kde_eldercare_bs, asp=1, ncol=3)
## OpenStreetMapData read with read_osm is static, so not usable in view mode. Please use tm_basemap or tm_tiles, with the provider name set to "OpenStreetMap.Mapnik"
## OpenStreetMapData read with read_osm is static, so not usable in view mode. Please use tm_basemap or tm_tiles, with the provider name set to "OpenStreetMap.Mapnik"
## OpenStreetMapData read with read_osm is static, so not usable in view mode. Please use tm_basemap or tm_tiles, with the provider name set to "OpenStreetMap.Mapnik"
tmap_arrange(kde_chas_jw, kde_infocomm_jw, kde_eldercare_jw, asp=1, ncol=3)
## OpenStreetMapData read with read_osm is static, so not usable in view mode. Please use tm_basemap or tm_tiles, with the provider name set to "OpenStreetMap.Mapnik"
## OpenStreetMapData read with read_osm is static, so not usable in view mode. Please use tm_basemap or tm_tiles, with the provider name set to "OpenStreetMap.Mapnik"
## OpenStreetMapData read with read_osm is static, so not usable in view mode. Please use tm_basemap or tm_tiles, with the provider name set to "OpenStreetMap.Mapnik"